function [p_lme, dp_lme, s_near]=ex33_wrapper_lme(x_n,x_s,options)
%function [p_lme, dp_lme, s_near]=ex33_wrapper_lme(x_n,x_s,options)
% INPUT:
%   x_n: node coordinates
%   x_s: sample points
%   options: parameters setting
%
% OUTPUT:
%   p_lme : are the shape functions values for each sample point
%   dp_lme: gradient for the shape functions at the sample points
%   s_near: samples adjacency, for each sample point this contains the indices of the
%           nearest nodes

if isfield(options, 'verb')
  verb = options.verb;
else
  verb = 1;
end
t=cputime;

if isfield(options, 'dim')
  dim = options.dim;
else
  error('Variable dim has not be defined in options');
end

%% samples adjacency structure is computed
beta_n  = options.beta_n;
range_n = options.range_n;
s_near  = SamplesAdjacency(x_n,x_s,range_n);


%% LME Shape functions and spatial derivatives up to second order computation
n_s = size(x_s,1);
p_lme = {zeros(n_s,1)};
dp_lme= {zeros(n_s,1)};
TolNR = options.TolNR;  %Newton-Raphson tolerance
for k=1:n_s
  nn_ids = s_near{k}; %nearest neighbour indices list
  n_near = length(nn_ids);
  x_near = x_n(nn_ids,:);

  %beta_n have the values of beta_a for each neighbor of xsample(i)
  beta_a = beta_n(nn_ids);

% shape functions, gradient and hessian
  [p_lme{k} dp_lme{k}] = ...
    shapef2_once(dim,TolNR,n_near,beta_a,x_near,x_s(k,:));
end

if verb==1
  format compact
  fprintf(1,'CPUTIME LME: %f\n', cputime-t);
end